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AFIT/GAI/ENY/03-3 


Abstract 

The  feasibility  of  using  a  microsatellite  to  accomplish  an  orbital  rendezvous  with  a 
noncooperative  target  is  being  evaluated.  This  study  focused  on  the  control  laws 
necessary  for  achieving  such  a  rendezvous.  The  relative  motions  of  the  microsatellite  and 
the  target  satellite  were  described  using  Hill’s  equations  and  two  different  controller 
methodologies  were  investigated.  An  impulsive  thrust  controller  based  on  the  Clohessey- 
Wiltshire  solution  was  found  to  use  little  fuel,  but  was  not  very  robust.  A  continuous 
thrust  controller  using  a  Linear  Quadratic  Regulator  was  found  to  be  more  robust,  but 
used  much  more  fuel.  As  a  final  solution,  a  hybrid  controller  was  evaluated  which  uses 
the  low  thrust  Clohessey-Wiltshire  approach  to  cover  most  of  the  necessary  distance,  and 
then  switches  to  the  Linear  Quadratic  Regulator  method  for  the  final  rendezvous  solution. 
Results  show  that  this  approach  achieves  rendezvous  with  a  reasonable  amount  of  control 
input. 


xi 


A  STUDY  OF  CONTROL  LAWS  FOR  MICROSATELLITE  RENDEZVOUS  WITH  A 

NONCOOPERATIVE  TARGET 


I.  Introduction 


Background 

The  United  States  derives  great  benefit  from  space-based  assets,  and  U.S.  dependence 
on  space-based  capabilities  will  only  continue  to  increase.  Other  countries  who 
recognize  the  advantage  conferred  on  the  U.S.  by  its  space  prowess  may  wish  to  develop 
the  means  to  neutralize  it. 

One  potential  method  would  be  to  use  a  parasitic  satellite.  A  parasitic  satellite  would 
rendezvous  with,  and  potentially  attach  to  a  target  satellite,  where  it  would  await  a 
command  from  the  ground  to  either  disrupt  satellite  operations  or  destroy  the  satellite. 

The  Chinese  claim  to  be  developing  microsatellites  to  perform  this  role.  (7)  The 
concept  of  operations  would  likely  begin  with  a  ground-based  orbit  determination  of  the 
target  satellite.  Then,  a  microsatellite  would  either  be  launched  or  released  from  a  mother 
ship  that  is  already  on  orbit.  Finally,  the  microsatellite  would  autonomously  achieve 
rendezvous  and  attach  to  the  target. 

Microsatellites  have  a  mass  of  100  kg  or  less,  and  are  less  expensive  to  build  and 
launch  than  larger  satellites.  However,  an  inherent  limitation  in  such  a  small  vehicle  is 
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fuel  capacity.  Therefore,  the  control  laws  used  to  achieve  rendezvous  should  minimize 


fuel  usage  to  the  maximum  extent  possible. 


Problem  Description 

In  an  effort  to  assess  the  threat  posed,  this  study  focuses  on  the  control  laws  necessary 
for  a  microsatellite  to  achieve  orbital  rendezvous  with  a  non-cooperative  target.  Other 
aspects  of  microsatellite  capabilities  are  being  studied  by  other  researchers.  Once  all  of 
the  various  investigations  have  been  completed,  it  is  expected  that  the  compilation  of 
results  will  indicate  the  overall  feasibility  of  the  proposed  system. 

For  this  project,  it  is  assumed  that  the  microsatellite  will  be  placed  into  an  orbit  similar 
to  that  of  the  target  satellite,  approximately  1000  km  behind  it  in  the  same  orbital  plane. 
The  microsatellite  then  performs  rendezvous  maneuvers  to  approach  the  target. 

For  this  project,  it  is  assumed  that  the  microsatellite  has  perfect  knowledge  of  the 
target’s  position  and  velocity  at  all  times.  In  reality,  the  microsatellite  would  likely  begin 
with  an  orbit  solution  derived  from  off-board  sensors.  As  the  microsatellite  approached 
the  target,  onboard  sensors  would  detect  the  target  satellite  and  an  updated  orbit  solution 
would  be  calculated.  This  would  allow  the  microsatellite  to  complete  the  rendezvous 
without  any  feedback  from  the  target  satellite.  Future  work  on  this  topic  should  include 
accounting  for  the  uncertainties  that  would  exist  in  reality. 

This  study  begins  with  satellite  dynamics;  then,  the  considered  control  methodologies 
are  discussed.  Finally,  the  controllers  are  employed  and  the  results  evaluated. 
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II.  Dynamics 


Coordinate  Frames 

Three  different  coordinate  frames  are  used  throughout  this  project.  The  first  is  the 
Earth  Centered  Inertial  (ECI)  frame,  which  is  depicted  in  Figure  1.  The  ECI  frame  is 
inertially  fixed  in  space  and  has  its  origin  at  the  center  of  the  Earth.  The  first  axis  in  the 
ECI  frame  points  toward  the  vernal  equinox,  the  second  axis  is  nonnal  to  the  first  in  the 
equatorial  plane,  in  a  direction  that  completes  the  frame  with  the  third  axis  pointing  out 
of  the  North  Pole. 
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Figure  1  Earth  Centered  Inertial  (ECI)  Coordinate  Frame 
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A  second  coordinate  frame  is  the  Perifocal  frame  (PQW).  The  PQW  frame  also  has 
its  origin  at  the  center  of  the  Earth,  but  the  first  two  axes  are  in  the  orbital  plane  of 
interest.  The  first  axis  points  toward  perigee  and  the  second  axis  is  perpendicular  to  the 
first  such  that  the  third  axis  points  in  the  direction  of  the  cross  product  of  the  satellite’s 
position  and  velocity  vectors.  This  frame  is  shown  in  Figure  2. 
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Figure  2  Perifocal  (PQW)  Coordinate  Frame 


The  third  coordinate  frame  used  in  this  project  is  Hill’s  (RTZ)  coordinate  frame. 

Hill’s  coordinate  frame  is  also  in  the  plane  of  the  orbit  of  interest,  but  it  includes  a 
reference  orbit  which  must  be  circular.  A  reference  origin  moves  about  the  circular  orbit 
with  mean  motion,  n  .  The  first  axis  points  in  the  radial  direction,  the  second  axis  points 
in  the  direction  of  the  instantaneous  velocity,  and  the  third  axis  points  in  the  out  of  plane 
direction  corresponding  to  the  cross  product  between  the  first  two  axes.  Hill’s 
Coordinate  Frame  is  illustrated  in  Figure  3. 
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Figure  3  Hill’s  (RTZ)  Coordinate  Frame 


It  is  important  to  be  able  to  transform  vectors  from  one  coordinate  frame  to  another. 
This  can  be  done  by  finding  a  transformation  matrix.  This  matrix  can  be  multiplied  by  a 
vector  in  one  coordinate  frame  to  find  its  components  in  another  frame.  Since  the 
Perifocal  frame  and  Hill’s  frame  are  both  in  the  same  orbital  plane,  the  transformation 
matrix  need  only  account  for  rotation  about  one  axis.  For  example,  to  transform  a 
velocity  vector  from  the  RTZ  frame  to  the  PQW  frame,  the  transformation  matrix,  S,  is 
shown  in  Equation  1,  where  v  is  the  angle  between  the  p  axis  and  the  er  axis. 


cosv 

sinv 

0 


-  sin  v  0]  [  v- 
cos  v  0  V.  =  S[vr&] 

0  lJU. 


(1) 
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This  case  is  illustrated  in  Figure  4,  where  w  and  ez  are  collinear,  out  of  the  page. 


q 


Figure  4  Transfonnation  From  RTZ  to  PQW  Coordinate  Frame 


Transforming  a  vector  from  the  PQW  frame  to  the  ECI  frame  can  include  up  to  three 
rotations,  involving  the  right  ascension  of  the  ascending  node,  Cl ,  the  argument  of 
perigee,  co ,  and  the  inclination,  /.  The  three  rotation  matrices  are  combined  to  give  a 
single  transfonnation  matrix,  R : 


R  = 


cos  Q  cos  co  -  sin  Q(cos  /) sin  co 
sin  Q  cos  co  +  cos  Q(cos  / )sin  co 
(sin/) sin  co 


-  cos  Cl  sin  co  -  sin  Q(cos z)cos  co 
-sinQsinry  +  cosQ(cos/)cos<y 
(sin/) cos  co 


sinQ(sinz') 

-cosQ(sinz') 

(cos/) 


(2) 


which  is  multiplied  by  a  vector  in  the  PQW  frame  to  find  its  components  in  the  ECI 
frame.  Equation  3  shows  an  example  for  transfonning  a  position  vector.  (9) 

[qJ=4Vv]  (3) 
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To  transform  a  vector  from  the  ECI  frame  to  the  PQW  frame,  the  inverse  of  the 


transformation  matrix  is  used,  such  as  in  Equation  4. 

\rPqw\  =  R%k]  (4) 

Similarly,  a  vector  can  be  transformed  from  the  PQW  frame  to  the  RTZ  frame  by  using 
the  inverse  of  the  S  transformation  matrix: 

[vr<k]  =  S~l[vpqw\  (5) 


Orbit  Characterization 

If  the  Earth  is  considered  to  be  inertially  fixed  in  space,  a  satellite  in  orbit  around  the 
Earth  has  six  degrees  of  freedom.  Therefore,  specifying  six  orbital  parameters 
completely  defines  the  orbit.  Two  sets  of  parameters  that  can  be  used  are  the  three 
dimensional  position  (f )  and  velocity  (v),  or  the  classical  orbital  elements  (COE),  which 
are  identified  in  Table  1.  (12) 


Classical  Orbital  Elements 

a 

Semi-Major  Axis 

e 

Eccentricity 

V 

True  Anomaly 

i 

Inclination 

Q 

Right  Ascension  of  the  Ascending  Node 

CO 

Argument  of  Perigee 

Table  1  Classical  Orbital  Elements 
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Both  sets  of  parameters  are  used  in  this  project,  and  it  is  convenient  to  be  able  to 
transfonn  from  one  set  to  the  other.  To  transform  from  the  position  and  velocity  vectors 
to  COEs,  the  first  step  is  to  calculate  the  angular  momentum,  H : 


H  =  r  xv 


(6) 


Next,  the  eccentricity  can  be  found: 


vxH  r 
ju  \r\ 

where  /u  is  the  Earth’s  gravitational  parameter  (ju  =  398601  hn  / s2  ).  (12)  The  semi¬ 
major  axis  can  be  calculated  with  Equation  8. 


(V) 


H 

2 

1- 

e 

1 

(8) 


The  inclination  is  a  function  of  the  angular  momentum  and  the  unit  vector,  k ,  where 
k  =  [0  0  1]: 


i  =  cos 

Next,  a  unit  vector,  n  ,  is  found  in  the  direction  of  the  ascending  node: 

„  k  x  H 


k*H 

H 


n  = 


kxH 


(9) 


(10) 


The  right  ascension  of  the  ascending  node  can  now  be  calculated  with  Equation  11,  and  a 
quadrant  check  perfonned: 
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Q  =  cos  1  (77  •  i ) 


(ID 


(n  •./)-  0  0  <  Q  <  ;r 

(77  •  j)  <  0  n  <  Cl  <271 


The  argument  of  perigee  is  found  from  Equation  12,  and  also  requires  a  quadrant  check: 


co  =  cos 


^  n  •  e  ^ 


V  |«|  J 


(12) 


(e  •  k)>  0  0  <  &>  <  ;r 

(e*k)<  0  n<(o<2n 


Finally,  the  true  anomaly  can  be  calculated,  along  with  a  quadrant  check: 

_,f  n»r^ 

V  =  cos  - 

l  \H  J 

(77  •  v)  <  0  0  <  v  <  ;r 

(77  •  v)  >  0  n  <v  <2n 


(13) 


To  transform  from  COEs  to  the  position  and  velocity  vectors,  the  magnitude  of  the 
position  and  velocity  are  found  using  Equation  14  and  Equation  15: 

?(l  ~|g|2) 


r  =  ■ 


1  +  lei  cos  v 


(14) 


pl=Jwi 

Next,  the  position  and  velocity  vectors  are  found  in  the  PQW  frame  as  follows: 


tSk 

1 _ 

cosv 

|  r 

sinv 

1 

O 

(15) 


(16) 
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(17) 


|v|(—  sin  v) 

|v|(|e|  +  cosv) 

0 

The  position  and  velocity  can  now  be  transfonned  into  the  ECI  frame  using  Equation  2 
and  Equation  3. 


Two-Body  Motion 

Satellite  motion  can  be  described  by  the  two-body  equations  of  motion,  where  the 
satellite  and  the  Earth  are  the  two  bodies  of  interest.  The  origin  in  Figure  5  is  inertially 

fixed,  Re  is  the  position  of  the  Earth  relative  to  the  origin,  Rs  is  the  position  of  the 
satellite  relative  to  the  origin,  and  r  is  the  position  of  the  satellite  relative  to  the  Earth. 
(9) 


Origin 


Earth  r  Satellite 

Figure  5  Two-Body  Motion 


r  can  be  found  from  RE  and  Rs  through  the  following  relation: 

r  =  Rs~RE  (18) 
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Differentiating  Equation  1  twice  yields: 


?  =  RS -Re 


Next,  the  forces  on  the  satellite  can  be  summed: 


=  msa  =  msRs 


The  only  force  acting  on  the  satellite  is  the  Earth’s  gravity: 


Fs= 


Gnuin, 


R s  RE 


where  ms  and  me  are  the  mass  of  the  satellite  and  the  mass  of  the  Earth,  respectively, 
and  G  is  the  universal  gravitational  constant.  The  unit  vector  is  defined  as: 


Rs-Re 
Rs  ~Re 


Equating  Equation  20  and  Equation  2 1 : 


Gm  c  m , 


__  V_»  l/l'  r/ft  n  /v 

mcR  c  = - -u 


\RS  Re 


Substituting  Equation  22  into  Equation  23: 


s  Gm^mF  -  - 

msRs  =  ~—  _  .-i  (Rs  —  Re) 


Rs-Re 


and  dividing  through  by  ms  yields  the  satellite’s  acceleration  with  respect  to  the  orgin: 


Rs=-t 


<Rs-Re) 


rs-re 


Similar  treatment  for  the  Earth  yields: 
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Rr  =■ 


Gnic 


Rs  Re 


ARs-Re) 


(26) 


Substituting  Equation  25  and  Equation  26  into  Equation  19: 

G(ms  +  mE) , 


r  =  —- 


Rs-Re 


■( Rs-Re ) 


(27) 


Assuming  ms «  me ,  let: 


G{ms  +  mE)  ~  GmE  =  // 


(28) 


Where  //  is  the  Earth’s  gravitational  parameter  (//  =  398601  kin' /.s2 ).  Finally: 


r  =  — 


jur 


(29) 


Perturbations 

Equation  29  describes  two-body  motion  without  any  other  effects  present,  such  as 
Earth  oblateness,  atmospheric  drag,  third  body  effects,  or  solar  wind.  Adding  a 
perturbation  tenn  to  Equation  29  yields: 

^  =  -7 ZT  +  ap  (30) 

H 

For  this  project,  the  J2  and  drag  perturbations  will  be  the  only  ones  considered: 

a  p  =  aJ2  +aD  (31) 
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Earth  Oblateness 


The  Earth  is  not  a  perfect  sphere;  rotational  motion  causes  it  to  bulge  about  the 
equator.  The  uneven  mass  distribution  produces  a  non-unifonn  gravity  field,  which 
causes  periodic  affects  to  a  satellite’s  orbital  elements.  Additionally,  the  Right  Ascension 
of  the  Ascending  Node,  Q ,  and  the  argument  of  perigee,  at ,  experience  secular  affects. 
These  effects  can  be  modeled  using  potential  theory.  (1) 

Expanding  r  into  its  constituent  x,  y,  z  components  in  the  /' ,  j  ,  k  directions: 


jux 


2  .  2  , 

\X  +  y  +  z  )2 


3  +aJ2i 


l  + 


■w 


2  ,  2  ,  2 

\x  +  y  +z  f 


T  +  aJj 


j  + 


fJZ 


2  ,  2  ,  2  \T 

\X  +y  +Z  Y 


T  +  aJ2k 


Next,  let  B  be  a  potential  function,  such  that 

VB  =  a , 


(32) 


(33) 


Therefore, 


d_ 

8x 


\ 

B 


-  jUX 

(x2  +y2  +z2y 


+  a 


J2i 


d_ 

dy 


\ 

+  B 

J 


W 


T  +  aJ2j 


/  2  .  2  ,  2  \  ; 

\X  +y  +z 


d_ 

dz 


^  +  B 


-/jz 


T  +  aJ2k 


(xz  +y~  +zz) 


2\2 


(34) 


(35) 


(36) 
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Comparing  Equations  34,  35,  and  36  with  32,  it  can  be  seen  that 


r  =  V 


f  \ 

±  +  B 

vH  j 


(37) 


where  B  is  an  infinite  series  that  models  the  Earth’s  non-homogeneous  mass  distribution: 

["/  \n  /  \n  1 

_  . .  co  D  n  D 

5  =  “fti2  ft  J»pn(s in^)+Z  ft  (c™cos^+‘s-sin^)/5™(sin^)  r  (38) 

T*  \V'\  \T*\ 

y  |  n= 2  \  h  |  /  /m=1  \  p  |  / 

and  <p  =  m2  ,  Re  is  the  Earth’s  radius  at  the  equator  ( Re  =  6378.135  km  ),  (f>  is  the 
latitude  measured  from  the  equator,  |r|  is  the  magnitude  of  the  satellite’s  position  vector, 
J „ ,  Cnm  and  Snm  are  the  zonal,  tesseral  and  sectorial  harmonic  coefficients,  respectively, 
and  P„  and  Pnm  are  Legendre  polynomials.  (2) 


k 


Figure  6  Geocentric  Spherical  Coordinate  Frame 
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The  Earth’s  oblateness  has  a  zonal  hannonic  effect,  symmetric  about  the  Earth’s  axis 
of  rotation  and  independent  of  longitude.  The  J2  contribution  is  much  larger  than  that  of 
other  harmonic  coefficients,  and  will  therefore  be  the  only  one  included  here.  Reducing 
Equation  38  and  considering  only  the  case  when  n  =  2  yields: 

_  ( R  y 

B  =  —^~  (39) 

H  UrU 

where  J2  =  0.0010826  and  the  second  Legendre  polynomial  is 

P2  (sin  <f)  =  ^  (3(sin  (ff  - 1)  (40) 

Substituting  Equation  40  into  Equation  39  and  rearranging: 

B  =  BJlRe  [|  ~  3(sin  (/)  f  ]  (41) 

2|r 

From  the  geometry  in  Figure  6: 


Substituting  Equation  42  into  Equation  41 : 


(42) 


(43) 
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taking  the  gradient: 


and  simplifying: 


(45) 


Drag 

Atmospheric  drag  can  alter  the  orbit  of  a  satellite  below  about  1000  km  in  altitude. 
The  drag  acceleration  is  added  to  the  two-body  equations  of  motion  as  a  second 
component  of  the  perturbation  term,  as  shown  in  Equation  3 1 .  The  drag  acceleration  is 
calculated  as  follows: 
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(46) 


aD=~TP 


V 


2  CdA 


m 


where  p  is  atmospheric  density,  V  is  the  magnitude  of  the  velocity  vector,  CD  is  the 


drag  coefficient,  A  is  the  satellite  cross-sectional  area  presented  to  the  atmosphere,  m  is 
the  mass  of  the  satellite,  and  iv  is  a  unit  vector  in  the  direction  of  the  velocity  vector. 


(1) 


The  density,  p ,  is  calculated  using  the  following  equation  and  a  standard  atmosphere: 


P  =  P„e 


(h-hQ) 

H 


(47) 


where  po  is  the  density  at  a  reference  altitude,  h  is  the  current  altitude,  hQ  is  the  reference 

altitude,  and  H  is  the  scale  height  for  the  reference  altitude.  (10) 

The  drag  coefficient,  Co  ,  usually  ranges  from  2.2  for  a  sphere  to  3.0  for  a  cylinder. 

( 1)  For  this  project,  the  drag  coefficient  for  the  target  satellite  was  set  at  2.2  and  for  the 
microsatellite  it  was  set  at  3.0.  This  represents  a  worst-case  scenario  for  the 
microsatellite’s  drag  coefficient. 

For  the  target  satellite,  the  cross-sectional  area  was  set  at  4.2  m  and  the  mass  was 
entered  as  725  kg.  These  parameters  correspond  to  Defense  Meterological  Satellite 
Program  (DMSP)  satellites.  They  were  selected  as  representative  of  an  operational 
satellite,  although  DMSP  satellites  orbit  at  higher  altitude  than  the  target  satellite  used  in 
this  project.  (5) 

A  mass  of  100  kg  and  a  cross-sectional  area  of  1.5  m  were  used  for  the  microsatellite. 
The  mass  was  selected  based  on  the  definition  of  a  microsatellite,  and  the  relatively  large 
cross-sectional  area  was  used  to  represent  worst-case  performance. 
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Relative  Motion 


For  rendezvous  maneuvers,  it  is  useful  to  know  the  positions  and  velocities  of  the 
microsatellite  and  the  target  relative  to  a  circular  reference  frame.  In  Figure  7,  O  is  the 
origin  centered  in  the  Earth  and  fixed  in  inertial  space.  O'  is  the  origin  of  a  reference 
frame  that  is  centered  on  the  instantaneous  location  of  a  point  moving  about  O  in  a 
circular  orbit  with  mean  motion,  n  .  The  unit  vectors  in  the  circular  reference  frame  are 
er,ee,  e.  in  the  radial,  in-track,  and  out  of  plane  directions,  respectively,  and  ro  is  the 

radius  of  the  circular  reference  orbit.  To  find  the  equations  for  relative  motion,  the  left 
side  of  Equation  29  is  found  by  taking  the  second  derivative  of  position,  and  the  right 
side  of  Equation  29  is  found  by  expanding  terms.  Then,  the  two  sides  are  equated.  (9) 
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Figure  7  Hill’s  Coordinates 


In  this  frame,  the  position  of  the  satellite  is: 

r  =  [(r0  +  8r)cos86^r  +  [(r0  +  (if)sin  86\en  +  [&]e,  (48) 

and  the  velocity  can  be  found  from: 

v=  —  (f)=  — (r)  +  (/7xr)  (49) 

dt  dt 

where  the  superscripts  i  and  o  correspond  to  the  inertial  and  circular  reference  frames, 
respectively,  and  the  mean  motion  of  the  circular  reference  frame  is: 

(50) 


Applying  Equation  49: 
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(51) 


v  =  \sdi'0  (-  sin  86)  +  Sr  (cos  86)  -  8r(  sin  86)86}  r 
+  ra  (cos  86)86  +  fir  (s  i  n  86)  +  8rSd( cos  86)}g 
+ 

+  (nez )  x  {[(ro  +  fir)cosfifi]fr  +  [(rc  +  fir)  sin  86\en  +  [&]<?_ } 


combining  terms: 

v  =  [fir  cos 86  -  [3d  +  n}ro  +  Sr)smSd}r  +  }S6  +  n\ra  +  Sr) cos 86  +  Sr  sin  86}0  +  [&]ez 

(52) 

The  position  (Equation  48)  and  velocity  (Equation  52)  in  the  relative  frame  can  now  be 
linearized  by  assuming  Sr,  86,  8z,  Sr,  86  and  8z  are  small.  Canceling  higher  order 
terms,  and  using  the  small  angle  approximations  cos 86  ~  1  and  sin  86  «  86  yields: 

r  ~  (r0  +  8r)e,.  +  (roS6)e0  +  (dz)e:  (53) 

v  «  (8r  -  nra86)er  +  \ra86  +  n(r0  +  8r)}g  +  (&)e,  (54) 

Now,  taking  the  derivative  of  velocity  yields  the  linearized  acceleration: 

«=  4-(v)=  ■^■(v)  +  (»xv)  (55) 

dt  at 

a  =  [Sr  -  nra8d)er  +  (; r086  +  n8r]eg  +  ( 8z)e ,  -  (nro86  +  n2 (ra  +  8r  f)er  +  (riSr  -  n2ra86)eg 

(56) 

combining  terms: 

a  =  r  =  \Sr  -  2nr0S6  -n2 (ra  +  8r )]f r  +  ( raS6  +  2 n8r  -  n2r0S6)eg  +  ( 8z)e:  (57) 

Equation  57  represents  the  left  side  of  Equation  29.  To  find  the  expression  for  the 
right  side  of  the  equation  in  the  relative  frame,  the  terms  are  expanded  using  the 
linearized  position  vector  (Equation  53).  In  the  denominator: 
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M3 


;  +  Sr)1  +  (roS8f  +  (&)2 1  =  (rQ2  +  2  ra8r  +  Sr 2  +  r2  SO1  +  Sz2 ) 


Sr  Sr'  2 
1  +  2 - 1 — —  +  SO  +  — 

ro  r0  rQ~ 


=  U 


V 


Sr  Sr'  2 
1  +  2 - 1 - — +  <%?  +  — 

ro  J 


2  (58) 


r  r 

O  '  O 


Therefore: 


( 

(  \ 

/ur 

\JL 

1  — 13  _ 

3 

r 

U  ) 

[fe  +  Sr)e,.  +  (ro<$0)e,  +  (Sz)ez  ] 


(59) 


-  Sr  Sr2  2 
1  +  2 - 1 - —  +  SO  +  — 


v 


O  / 


A  binomial  expansion  can  be  used  for  the  denominator:  (3) 
1 


,  Sr  Sr'  2 
1  +  2  —  +  —  +  S0-+  — 


v 


r  r 

O  '  o 


3 

f  Sr  Sr2  2 

1 - 

2  +  .  +  SO2  +  _ 

2 

N. 

k: 

'o  J 


15 

H - 

8 


Sr  Sr'  2  & 

2 +  — - — h  SO  + — ; 


2  b 


V  ro 


+  Higher  Order  Terms 


O  J 


(60) 

It  can  be  seen  that  the  second  order  and  higher  terms  will  not  be  linear,  and  will  therefore 
be  cancelled.  This  yields: 

1 


f 


Sr  Sr 2  2 


1-3— 

r. 


(61) 


1  +  2 - 1 — - — I-  SO  +  0 

V  ro  ro~  Vo  J 


Recognizing  the  magnitude  of  the  mean  motion  from  Equation  32  and  substituting 
Equation  61  into  Equation  59: 

Sr 


=  ~n  2 1 fc  +  Sr)er  +  (roS6)ee  +  (&)ez  ] 


1-3- 


(62) 
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Multiplying  and  again  eliminating  higher  terms  produces  the  final  result  for  the  right  side 
of  Equation  29: 


—  7TT  =  -»2  [(>;  -  2 Sr)er  +  {roS6)e0  +  {Sz)e: ] 


Equation  57  and  Equation  63  can  now  be  equated: 


Sr  -  2nroS6  -  n 2  {ro  +  Sr)  =  -n2  (ro  -  2 Sr) 
roS0  +  2n Sr  -  n2roS0  =  -n2roS0 
Sz  =  -n2Sz 


Solving  Equation  64  yields: 


and: 


Sr(t) 

3  n  sin  i//  0 

0 

<5>(0) 

roS6(t) 

= 

6n(cos^-l)  0 

0 

’■Mo) 

&{t) 

0  0 

-  n  sin  yr 

5z(  0) 

cos^ 

2sin^ 

0 

Sr{0) 

+ 

-  2  sin  yr 

4cos^  -3 

0 

•■Mo) 

0 

0 

cos  y/ 

&(  o) 

Sr(t) 

4-3cos^  0 

0 

(>/•(() ) 

roS0{t) 

= 

6(sin^-^)  1 

0 

rM<>) 

Sz{t) 

0  0 

cos  y/ 

&(o) 

smy/ 


2-2cos^ 


2cos^"-2  4sin^-3^ 


n 

0 


n 

0 


0 

0 

smy/ 

n 


fif(O)  ' 

rM  0) 

Sz(  0) 


(63) 


(64) 


(65) 


(66) 


where  y/  =  nt .  In  compact  notation: 

[«(»)]  =  0,r[^(0)]+®„[«(0)]  (67) 
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and: 


\Sr  (t )]  =  <t>  rr  \dr  (0 )]  +  <t>  n,  [<5v  (0 )]  (68) 

These  equations  are  only  valid  when  dr  and  dz  are  small,  although  rod0  can  be  large. 
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III.  Control 


Both  of  the  control  methodologies  evaluated  in  this  thesis  are  based  on  linear 
equations  of  motion,  yet  the  orbits  are  propagated  using  non-linear  equations.  The  linear 
equations  are  only  valid  near  a  circular  orbit,  and  normally  are  used  with  a  leader  satellite 
in  a  perfectly  circular  orbit.  The  follower  satellite  remains  “close”  to  the  leader,  and  its 
position  and  velocity  are  found  relative  to  the  leader  satellite.  The  equations  are 
linearized  about  a  circular  lead  orbit,  so  the  lead  satellite’s  orbit  must  be  circular  for 
proper  application  of  the  equations.  For  this  project,  it  was  decided  to  allow  the  target 
satellite  to  have  some  eccentricity  in  its  nominal  orbit.  The  microsatellite  will  have 
eccentricity  introduced  into  its  orbit  as  the  control  solutions  are  applied.  To  allow  for 
both  satellites  to  have  non-zero  eccentricity,  and  to  still  comply  with  the  requirement  for 
the  leader  satellite  to  be  in  a  perfectly  circular  orbit,  it  was  decided  to  use  a  “virtual” 
leader.  The  virtual  leader  orbit  is  constructed  to  have  an  initial  position  vector  that  is 
collinear  with  the  target  satellite,  to  be  in  a  coplanar  orbit  with  the  target  satellite,  to  have 
the  same  semi-major  axis  length  and  therefore  the  same  period  as  the  target,  but  to  have 
zero  eccentricity.  The  virtual  reference  orbit  will  therefore  remain  close  to  the  target 
satellite  throughout  its  orbit.  The  position  and  velocity  of  the  target  and  the  microsatellite 
are  found  relative  to  the  virtual  circular  reference  orbit.  Instead  of  effecting  rendezvous 
to  the  origin  of  the  relative  circular  reference  frame,  the  microsatellite  is  able  to  pursue 
the  target  satellite  which  is  moving  with  respect  to  the  reference  frame.  (12) 


24 


For  both  of  the  control  methods  considered,  it  is  necessary  to  find  the  position  and 
velocity  of  a  satellite  relative  to  the  circular  reference  frame.  First,  the  position  and 
velocity  are  transfonned  into  the  circular  reference  orbit’s  PQW  coordinate  frame.  Then 
the  true  anomaly  of  the  reference  orbit  can  be  found  from  the  p  and  q  components  of  the 
reference’s  position  vector.  Using  MATLAB’s  atan2  function  ensures  quadrant 
accuracy: 


v  =  atan2 


\r P  ) ref 


(69) 


The  true  anomaly  can  then  be  used  to  set  up  a  transformation  matrix,  S,  which  is  used  to 
transfonn  the  satellite’s  position  from  the  circular  reference’s  PQW  coordinate  frame  to 
its  RTZ  coordinate  frame. 


(T'&  ) satellite  T pqw  ) 


P4W  '  satellite 


(70) 


The  atan2  function  can  now  be  used  to  find  86  without  quadrant  ambiguity: 


86  -atan2 


rr\ 

' e 


(71) 


\^r  J satellite 

The  position  vector  of  the  satellite  relative  to  the  circular  reference  orbit  is  calculated  as 
follows: 


8r  = 


(72) 


To  find  the  relative  velocity  vector,  it  is  useful  to  find  the  angle  between  the  reference 
orbit’s p  axis  and  the  satellite’s  position  vector  in  the pq  plane. 
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y  =  atan2 


(73) 


\rp  J 


satellite 


This  angle  can  be  used  to  build  a  transformation  matrix,  T,  which  will  be  used  to  find  the 
radial,  tangential,  and  out-of-plane  components  of  the  velocity  vector: 


(^123  ) 


satellite 


' ^  satellite 


(74) 


Next,  the  magnitude  of  the  satellite’s  angular  velocity  can  be  found: 

ife)* 


e 


/ satellite 


(75) 


Finally,  the  relative  velocity  can  be  found: 


(f) 


satellite 


e 


-n 


fe). 


satellite 


(76) 


where  n  is  the  mean  motion  of  the  circular  reference  orbit. 


Clohessey-Wiltshire 

The  first  control  methodology  considered  for  this  project  utilizes  the  Clohessey- 
Wiltshire  (CW)  equations.  This  technique  produces  an  impulsive  control  thrust  which 
initiates  the  rendezvous.  A  second  impulsive  thrust  is  applied  once  rendezvous  is 
complete,  to  null  out  the  relative  velocity  between  the  microsatellite  and  the  target 
satellite. 

Equation  68  is  rearranged  to  solve  for  the  relative  velocity  needed  at  time  t  =  0  in 
order  to  achieve  rendezvous  at  a  specified  time,  t  =  T  : 
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M>)]=  fc./PHr)]-®,,  [<?(0*  (77) 


To  calculate  the  Av  required  at  time  t  =  0 ,  the  necessary  inputs  are  the  microsatellite’s 
starting  position  relative  to  the  circular  reference,  the  desired  rendezvous  time,  the 
desired  position  of  the  microsatellite  relative  to  the  circular  reference  at  the  rendezvous 
time,  and  the  mean  motion  of  the  circular  reference  frame.  (12) 

The  desired  relative  position  for  the  microsatellite  at  the  rendezvous  time  is  the  same 
as  the  position  of  the  target  satellite  at  the  rendezvous  time.  The  target’s  initial  position 
and  velocity  relative  to  the  circular  reference  frame  can  be  found  as  described  above,  and 
then  Equation  68  can  be  employed  to  find  the  target’s  position  relative  to  the  circular 
reference  at  the  specified  rendezvous  time. 

The  microsatellite’s  initial  position  relative  to  the  circular  reference  frame  is  also 
found  using  the  procedure  described  earlier  in  this  chapter.  The  mean  motion  of  the 
circular  reference  frame  is  calculated  using  Equation  50.  Once  all  of  the  input  values  are 
entered  into  Equation  77,  the  necessary  relative  velocity  at  time  t  =  0  can  be  calculated. 
Equation  52  is  used  to  transform  the  velocity  from  the  relative  frame  to  an  inertial  frame. 
It  is  then  transformed  into  the  ECI  frame  and  applied  to  the  microsatellite,  which  is 
propagated  using  the  non-linear  equations  of  motion.  To  calculate  the  necessary  Av ,  the 
microsatellite’s  initial  velocity  is  subtracted  from  the  velocity  calculated  for  rendezvous. 


Av  =  (v...  )  —  (v... ) 

V  lJk  ' calculated  V  lJk  ' 


Vk  '  initial 


(78) 


At  the  specified  rendezvous  time,  a  second  impulsive  thrust  can  be  applied  to  null  the 
relative  velocity  between  the  microsatellite  and  the  target.  This  second  Av  is  added  to 
the  first  to  find  the  total  Av  for  the  CW  rendezvous. 
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The  Clohessey  Wiltshire  equations  can  be  improved  by  adding  the  J2  perturbation  to 
the  calculation,  which  is  discussed  in  the  next  section.  Although  drag  was  added  to  the 
equations  of  motion  in  the  previous  chapter,  the  drag  perturbation  is  not  accounted  for  in 
the  Clohessey  Wiltshire  equations  for  this  project. 


Clohessey  Wiltshire  with  J2 


McLaughlin,  et  al,  discuss  addition  of  the  J2  perturbation  term  to  Hill’s  equations.  (4) 
The  primary  change  is  substitution  of  core  in  place  of  n  in  the  in-plane  components 

{dr,  ro80 )  and  0):  in  place  of  n  in  the  out  of  plane  component  (&)  in  the  ®  matrices  in 
Equation  66,  where: 


G>rg  =n  +  M 

(79) 

a>_  =  n  +  M  +  d> 

(80) 

and: 


M  = 


fl>=3  nRE2j2(4-5sm2i) 

4  p2  V  ' 

where  M  is  the  effect  of  J2  on  the  mean  anomaly,  cb  is  the  effect  of 
of  perigee,  n  is  the  mean  motion,  RE  is  the  mean  radius  of  the  Earth 
the  perturbation  term,  e  is  the  eccentricity,  i  is  the  inclination  and  p 
rectum.  (4) 


(81) 

(82) 

J 2  on  the  argument 
at  the  equator,  J2  is 
is  the  semilatus 
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Further  terms  are  discussed  by  McLaughlin,  et  al,  to  account  for  the  J2  effects 
resulting  from  the  target  and  chase  satellites  having  different  inclinations.  However,  the 
inclinations  are  assumed  to  be  the  same  for  this  project,  so  the  extra  terms  are  not 
discussed  here. 

After  cor0  and  co,  are  substituted  into  Equation  66,  by  replacing  i//  =  nt  with 
Vrg  =a>rd 1  'n  the  in-plane  components  of  the  <f>  matrices,  and  replacing  i//  =  nt  with 
y/z  =  a>:t  in  the  out  of  plane  components  of  the  ®  matrices,  then  Equation  66  can  be 
rearranged  into  the  form  of  Equation  77  to  give  a  solution  for  the  necessary  velocity  as 
discussed  in  the  previous  section. 


Linear  Quadratic  Regulator 


The  second  control  methodology  considered  for  this  project  was  a  Linear  Quadratic 
Regulator  (LQR).  Again  utilizing  the  relative  reference  frame  (Figure  7)  and  the  relative 
equations  of  motion  (Equation  64),  a  vector  comprising  both  the  relative  position  and 
velocity  can  be  defined:  (2) 


Sr 

Sr 

rose 

roS0 

Sz 

Sz 


(83) 


and  its  derivative  is: 
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X  = 


(84) 


dr 
Sr 

rose 
rose 

Sz 
Sz 

The  relative  equations  of  motion  can  be  placed  in  state  equation  fonn: 

x  =  Ax  +  Bu  (85) 

where  u  is  a  vector  of  control  inputs.  Equations  64,  83  and  84  can  now  be  used  to  write 
Equation  85  as: 
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(86) 


A  Linear  Quadratic  Regulator  obtains  the  optimal  gain  matrix,  K  for  the  control 
vector:  (6) 


u  =  -Kx 


(87) 


by  minimizing  the  performance  index: 

oo 

J  =  j  {x' Qx  +  u'Ru)dt 

0 


The  algebraic  Riccati  equation  is  solved  to  find  S: 

SA  +  A'S-SBR1B'S  +  Q  =  0 


(88) 


(89) 
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where  Q  is  the  state  weighting  matrix  and  R  is  the  control  weighting  matrix.  Higher 


values  in  the  Q  matrix  speed  movement  toward  the  desired  state,  and  higher  values  in  the 
R  matrix  reduce  control  usage.  For  this  project,  Q  and  R  will  have  the  following  fonns: 

q  0  0  0  0  0 

0  q  0  0  0  0 

0  0  q  0  0  0 

0  0  0  q  0  0 

0  0  0  0  q  0 

0  0  0  0  0  q 

r  0  0 

R=  0  r  0  (91) 

0  0  r 

Finally,  the  optimal  gain  matrix  is  calculated: 

K  =  R'B'S  (92) 

For  this  project,  the  optimal  gain  matrix  is  found  by  entering  the  A,  B,  Q  and  R 
matrices  into  MATLAB’s  LQR  function.  Once  calculated  during  a  run,  the  gain  matrix 
never  needs  to  be  recalculated. 


For  this  project,  the  microsatellite  is  chasing  the  target  satellite  rather  than  the  circular 
reference.  Therefore,  the  control  input  is  based  on  the  difference  between  the 
microsatellite’s  state  vector  and  the  target’s  state  vector  (i.e.,  the  error  state): 


rSO . 


(93) 
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The  following  steps  are  repeated  during  propagation.  First,  the  position  and  velocity 
of  the  target  and  microsatellite  are  found  relative  to  the  circular  reference.  Then  the 
control  input  is  calculated  using  Equation  93.  The  control  input,  u  ,  is  actually  a  specific 
thrust,  so  the  calculated  value  is  already  in  an  inertial  frame.  It  is  then  transformed  into 
the  ECI  frame  and  added  to  the  two-body  equations  of  motion: 

f  =  -zr  +  3p+u  (94) 

|r| 

The  target,  microsatellite  and  reference  orbits  are  propagated  one  time  step  using 
numerical  integration.  The  Av  for  each  step  is  calculated  using  Euler  integration;  the 
magnitude  of  the  control  acceleration  is  multiplied  by  the  size  of  the  time  step  and 
accumulated  over  the  run. 

The  next  iteration  begins  with  finding  the  position  and  velocity  of  the  target  and 
microsatellite  relative  to  the  circular  reference  frame.  The  process  is  repeated  until 
rendezvous  is  achieved.  This  algorithm  is  depicted  in  Appendix  B. 
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IV.  Results 


The  results  in  this  chapter  were  obtained  by  starting  with  a  target  orbit  having  the 
initial  orbital  elements  found  in  Table  2,  unless  otherwise  specified. 


Orbital  Elements 

Values 

Units 

a 

Semi-Major  Axis 

6772.888912204840 

kilometers 

e 

Eccentricity 

0.00098877135498 

dimensionless 

V 

True  Anomaly 

0.79736386485827 

radians 

i 

Inclination 

0.90757990078380 

radians 

Q 

Right  Ascension  of  the  Ascending  Node 

1.51843760980691 

radians 

CO 

Argument  of  Perigee 

5.59054044657763 

radians 

Table  2  Initial  Orbital  Elements  for  the  Target  Orbit 


The  starting  position  of  the  microsatellite  is  dr  =  [0  - 1000  0]  unless  otherwise 
specified. 


Simulator  Capability 

Orbits  based  on  the  non-linear  equations  of  motion  are  propagated  for  this  project 
using  numerical  integration.  The  explicit  Runge-Kutta  4,5  (Dormand,  Prince  pair) 
method  is  employed  in  MATLAB  using  the  ODE45  command.  For  this  project,  the 
relative  tolerance  has  been  set  to  2.2205e-14  and  the  absolute  tolerance  has  been  set  to 
le-200.  A  tighter  relative  tolerance  was  found  to  be  beyond  the  capability  of  the  system. 
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With  these  tolerances  set,  a  circular  orbit  was  propagated  ten  revolutions  using  the  non¬ 
linear  equations  of  motion;  the  same  orbit  was  propagated  ten  revolutions  using  the  linear 
equations  of  motion.  Since  the  circular  orbit  has  no  offsets  from  the  identical  circular 
reference  orbit,  the  criteria  for  linearization  are  perfectly  satisfied.  However,  when  the 
relative  position  between  the  linear  and  non-linear  solutions  is  compared,  a  small  error  is 
observed.  This  error  is  attributable  to  the  precision  of  the  simulator,  and  is  depicted  in 
Figure  8: 


..-11 

x  10 


Figure  8  Simulator  Error  -  Distance  Between  Linear  and  Non-Linear  Propagations 
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Linear  Versus  Non-Linear  Equations  of  Motion 


As  discussed  in  previous  chapters,  control  inputs  for  this  study  are  calculated  using 
methods  based  on  linear  equations  of  motion.  The  orbits  are  then  propagated  using  non¬ 
linear  equations  of  motion.  Before  discussing  the  results  of  this  project,  it  is  enlightening 
to  further  explore  the  influence  of  linearization. 

For  the  linearized  equations  to  remain  valid,  only  small  Sr  and  Sz  are  allowed; 
however,  roS6  may  be  large.  Adding  eccentricity  to  an  orbit  causes  it  to  have  a  greater 

variation  in  Sr  compared  to  the  circular  reference  orbit.  For  an  orbit  with  the  orbital 
elements  listed  in  Table  2,  there  is  a  clear  difference  in  the  orbits  propagated  using  the 
linear  and  non-linear  equations  of  motion.  The  results  are  seen  in  Figure  9  which  shows, 
in  the  Sr,  roS0  plane,  the  relative  distance  between  the  location  of  the  satellite  calculated 

using  the  linear  equations  of  motion  and  the  location  of  the  satellite  calculated  using  the 
non-linear  equations  of  motion. 
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Figure  9  Distance  Between  Linear  and  Non-Linear  Propagations  With  Eccentricity 

It  can  be  seen  that  there  is  a  periodic  difference  in  the  radial  direction  during  each 
orbit  and  a  secular  difference  in  the  in-track  direction.  The  peak  Sr  offset  in  each  orbit  is 
approximately  24  meters  and  the  roS6  accumulation  is  about  120  meters  per  orbit. 

Repeating  these  calculations  with  the  eccentricity  changed  to  0.00313260666343051 
produces  a  graph  of  similar  shape,  but  the  magnitudes  of  the  Sr  and  roS6  offsets  are 

approximately  240  meters  and  1 .2  kilometers  per  orbit,  respectively.  This  corresponds  to 
the  orbit  of  the  microsatellite  after  applying  the  delta-V  calculated  using  the  Clohessey- 
Wiltshire  equation.  Clohessey-Wiltshire  results  will  be  discussed  in  greater  detail  later. 

To  further  illustrate  the  differences  between  the  linear  and  non-linear  equations  of 
motion,  Figure  10  shows  the  difference  between  the  two  as  measured  in  an  inertial  frame 
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Total  (km 


and  Figure  1 1  shows  the  same  difference  as  measured  in  the  relative  frame.  Both  figures 
represent  the  difference  for  the  orbit  defined  in  Table  2  and  have  been  propagated  over 
ten  orbits. 


i 
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Figure  10  Difference  in  Linear  and  Non-Linear  Propagation  in  Inertial  Frame 
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The  secular  effect  in  the  ro86  direction  is  clearly  seen  in  the  second  graph  of 


Figure  1 1 . 
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Figure  1 1  Difference  in  Linear  and  Non-Linear  Propagation  in  Relative  Frame 


Clohessey-Wiltshire  Controller 

The  Clohessey-Wiltshire  equations  have  peculiar  behavior,  and  careful  selection  of  the 
desired  rendezvous  time  is  important  for  their  performance.  The  top  graph  in  Figure  12 
was  generated  by  varying  the  rendezvous  time  up  to  6  six  orbits.  The  reference  orbit’s 
semi-major  axis  length  was  used  to  calculate  the  necessary  mean  motion,  n  ,  the  starting 
relative  position  was  set  at  dr  =  [0  - 1000  0],  and  the  necessary  relative  velocity,  Sv  , 
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for  a  rendezvous  at  the  origin,  Sr  =  [0  0  0],  was  calculated.  The  magnitude  of  the 
calculated  velocity,  |<5v| ,  is  shown. 


Time  (orbits) 


Figure  12  Magnitude  of  the  Relative  Velocity  for  Rendezvous  to  the  Origin 


The  graph  reveals  that  some  selections  of  rendezvous  time  are  very  costly.  Between 
each  of  the  extreme  occurrences  are  a  range  of  more  reasonable  values.  The  lower  graph 
in  Figure  12  shows  a  closer  view  of  the  curve  in  the  range  of  rendezvous  times  from  3.6 
orbits  to  4.4  orbits.  Similar  troughs  are  seen  between  the  other  extreme  values.  As  the 
selected  rendezvous  time  increases,  the  troughs  reach  lower  \Sv\  values.  It  appears  that 

for  rendezvous  to  the  origin,  the  best  selection  of  rendezvous  times  are  integer  numbers 
of  orbits. 
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Since  the  target  in  this  study  has  non-zero  eccentricity,  it  is  unlikely  to  be  at  the 
reference  orbit’s  origin  at  rendezvous  time.  To  represent  the  potential  offset  from  the 
origin,  the  runs  from  Figure  12  were  repeated,  but  the  target  rendezvous  position  was  set 
to  dr  =  [lO  10  O] .  The  results  are  depicted  in  Figure  13.  The  closer  view  in  the  lower 
graph  reveals  the  effect  of  an  offset  target  is  to  add  additional  peak  values  in  the  troughs. 
These  new  peaks  appear  at  integer  numbers  of  orbits.  Therefore,  the  best  rendezvous 
times  for  this  case  are  within  each  trough,  close  to  an  integer  number  of  orbits  but  not  an 
exact  integer  number  of  orbits. 


0  I _ l _ l _ l _ _1 _ I _ I _ : _ I _ : _ I _ I _ I 
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Figure  13  Magnitude  of  the  Relative  Velocity  for  Rendezvous  to  an  Offset  Location 


Figure  14  shows  the  AV  curve  from  330  minutes  to  400  minutes,  for  the  rendezvous 
scenario  considered  in  this  chapter.  Since  the  target  satellite’s  orbital  period  is  92.45 
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minutes,  exactly  four  orbits  would  take  369.8  minutes.  A  rendezvous  time  of  368 
minutes,  which  is  slightly  less  than  four  orbits,  has  been  selected.  This  rendezvous  time 
requires  less  AV  than  is  found  in  troughs  associated  with  a  lower  number  of  orbits,  and 
within  the  trough  depicted  in  Figure  14,  it  corresponds  to  the  minimum  AV  immediately 
before  the  spike.  Using  a  rendezvous  time  from  a  trough  associated  with  a  higher  number 
of  orbits  would  reduce  the  required  AV ,  but  the  impact  of  propagating  a  linear  solution 
using  the  non-linear  equations  of  motion  is  secular,  resulting  in  a  larger  final  distance 
when  a  longer  rendezvous  time  is  used. 
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Figure  14  Delta-V  Required  to  Achieve  Rendezvous  for  the  Considered  Scenario 


Using  the  specified  rendezvous  time  of  368  minutes  and  including  no  perturbations, 
produces  the  results  in  Table  3.  The  initial  AV  is  required  to  begin  the  rendezvous 
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maneuver  and  the  second  AV  is  performed  at  the  rendezvous  time  to  null  out  the  relative 


velocity  between  the  microsatellite  and  the  target  satellite.  Since  rendezvous  to  the  target 
satellite  has  not  been  completely  accomplished,  further  maneuvers  will  be  necessary. 
Therefore,  matching  the  target’s  velocity  exactly  is  not  actually  desirable,  but  it  is 
accomplished  here  to  complete  the  results  from  a  single  rendezvous  maneuver. 
Subsequent  maneuvers  to  close  the  final  distance  using  the  Clohessey-Wiltshire 
controller  were  ineffective  at  completing  the  rendezvous  to  an  acceptable  distance, 
particularly  if  there  was  a  significant  Sr  offset.  This  can  be  readily  seen  in  the  examples 
illustrated  in  Tables  15  -  17  in  Appendix  A. 


CW  Rendezvous  -  No  Perturbations 

AFj  (km/sec) 

0.016539 

AV2  (km/sec) 

0.015269 

mai  (km/sec) 

0.031809 

Final  Distance  (km)  at  368  Minutes 

5.029880 

Table  3  CW  Rendezvous  -  No  Perturbations 


Figure  15  shows  how  the  relative  distance  between  the  microsatellite  and  target 
satellite  decreases  over  time.  The  distances  in  the  figure  were  calculated  in  the  relative 
reference  frame  and  propagated  with  the  linear  equations  of  motion. 
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Figure  15  Relative  Distance  Between  the  Target  and  Microsatellite  in  Relative  Frame 


The  same  Delta-V  solution  was  applied  and  propagated  in  the  inertial  reference  frame 
using  the  non-linear  equations,  and  the  results  of  both  propagations  are  depicted  in 
Figure  16  along  with  the  difference  between  the  two.  Both  graphs  have  similar  shape, 
but  the  difference  between  them  is  seen  to  be  both  periodic  and  secular  with  offsets  as 
large  as  23  kilometers  during  the  maneuver. 
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Figure  16  Linear  and  Non-Linear  Propagation  of  Cloehessey-Wiltshire  Solution 


In  Figure  17,  the  final  rendezvous  behavior  is  more  closely  examined.  It  can  be  seen 
that  the  linear  propagator  goes  to  a  relative  distance  of  zero  precisely  at  368  minutes,  as 
expected.  However,  propagation  using  the  non-linear  equations  has  the  effect  of 
offsetting  the  time  of  closest  approach,  to  a  distance  of  0.793406  kilometers  at  362 
minutes.  At  the  specified  rendezvous  time  of  368  minutes,  the  microsatellite  is  5.029880 
kilometers  from  the  target  satellite  and  is  already  moving  away.  This  figure  further 
illustrates  the  effect  of  using  linear  control  equations  to  calculate  maneuvers  which  are 
then  applied  using  non-linear  equations. 
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Figure  17  Closer  View  of  Final  Rendezvous  Behavior 

Adding  the  J2  perturbation  to  both  the  propagator  and  the  Clohessey-Wiltshire 
equations  produces  the  results  in  Table  4.  This  rendezvous  is  depicted  in  Figure  18, 
which  is  rendered  in  the  relative  coordinate  frame.  Since  the  target  has  a  non-zero 
inclination,  the  J2  perturbation  causes  out  of  plane  effects  during  the  rendezvous  which 
can  be  seen  as  divergence  in  the  dz  plot. 


CW  Rendezvous  -  J2  Perturbation 

AFj  (km/sec) 

0.016578 

AV2  (km/sec) 

0.015917 

mai  (km/sec) 

0.032496 

Final  Distance  (km)  at  368  Minutes 

3.518942 

Table  4  Clohessey-Wiltshire  Rendezvous  With  J2  Perturbation 
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Finally,  the  drag  perturbation  is  added  to  the  propagator.  This  produces  the  results 


listed  in  Table  5. 


CW  Rendezvous  -  J2  and  Drag  Perl 

turbations 

AVl  (km/sec) 

0.016578 

AV2  (km/sec) 

0.019170 

AVTo,ai  (km/sec) 

0.035748 

Final  Distance  (km)  at  368  Minutes 

3.321252 

Table  5  Clohessey-Wiltshire  Rendezvous  With  /?  and  Drag  Perturbations 


The  position  of  the  microsatellite  relative  to  the  target  satellite  is  depicted  in 
Figure  20,  for  the  Clohessey-Wiltshire  solution  with  J2  and  drag  perturbations  included. 


roSe  (km) 

Figure  20  Clohessey-Wiltshire  Rendezvous  with  Perturbations  in  Relative  Frame 
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Figure  2 1  Closer  View  of  Clohessey-Wiltshire  Rendezvous  with  J2  and  Drag 

A  review  of  the  distances  for  the  Clohessey-Wiltshire  cases  discussed  above  reveals 
that  the  microsatellite  is  closer  to  the  target  satellite  at  the  specified  rendezvous  time 
when  the  J2  perturbation  is  added  and  even  closer  when  drag  is  added.  A  closer  study  of 
the  results  reveals  that  rendezvous  time  is  shifted  from  the  specified  time  in  each  of  these 
cases,  and  the  closest  approach  does  occur  in  the  case  when  no  perturbations  are  present. 
Adding  perturbations  causes  the  closest  approach  distance  to  increase. 

Figure  22  shows  that  using  the  linear  equations  in  the  case  when  no  perturbations  are 
present  results  in  rendezvous  to  zero  distance  at  the  specified  time.  Propagating  the  same 
solution  using  the  non-linear  equations  causes  the  rendezvous  to  occur  6  minutes  early,  at 
362  minutes.  Adding  the  J2  perturbation  then  delays  the  rendezvous  by  3  minutes  from 


48 


the  non-perturbed  case,  to  365  minutes.  Addition  of  drag  further  delays  the  rendezvous 
to  370  minutes. 


CW  Rendezvous  Distances 

Case 

Distance  at  368 
Minutes  (km) 

Time  of  Closest 
Approach  (min) 

Distance  of  Closest 
Approach  (km) 

No  Perturbations 

5.029880 

362 

0.793406 

Ji  Perturbation 

3.518942 

365 

1.499542 

Ji  and  Drag 
Perturbations 

3.321252 

370 

2.398979 

Table  6  Clohessey-Wiltshire  Rendezvous  Distances 


Time  (minutes) 


Figure  22  Clohessey-Wiltshire  Rendezvous  Results 
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Linear  Quadratic  Regulator  Controller 


With  the  Clohessey-Wiltshire  controller,  AV  was  influenced  by  adjusting  the 
rendezvous  time.  With  the  Linear  Quadratic  Regulator  (LQR)  controller,  AV  is 
influenced  by  adjusting  the  values  in  the  state  weighting  and  control  weighting  matrices. 
Decreasing  control  usage  generally  saves  fuel,  but  rendezvous  may  be  delayed  or  not 
achieved. 

Since  rendezvous  time  is  not  fixed  for  this  method,  the  controller  is  allowed  to  operate 
until  the  criteria  that  define  a  successful  rendezvous  are  satisfied.  For  this  project,  two 
criteria  were  selected  that  must  be  concurrently  satisfied  for  a  successful  rendezvous: 

1 .  Achieve  a  relative  distance  of  1  meter  or  less 

2.  Achieve  a  relative  velocity  of  1  cm/sec  or  less. 

For  the  case  with  no  perturbations,  the  results  are  captured  in  Table  7.  The  state 
weighting  matrix  values  are  set  to  1  and  the  control  weighting  matrix  values  are  set  to 
1  e  1 3  for  this  case. 


LQR  Rendezvous  -  No  Perturbations 

^VTotal  (km/sec) 

0.306821 

Rendezvous  Time  (min) 

510 

Table  7  Linear  Quadratic  Regulator  Rendezvous  With  No  Perturbations 
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The  relative  distance  between  the  microsatellite  and  the  target  satellite  during  the 
rendezvous  is  shown  in  Figure  23.  The  figure  shows  that  the  majority  of  the  distance  is 
reduced  early  in  the  pursuit,  and  most  of  the  time  is  spent  eliminating  a  relatively  small 
final  distance. 
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Figure  23  Relative  Distance  During  LQR  Rendezvous  With  No  Perturbations 


To  compare  similar  rendezvous  maneuvers  using  both  of  the  controller  methodologies 
and  no  perturbations,  the  rendezvous  distance  and  velocity  criteria  for  a  Linear  Quadratic 
Regulator  rendezvous  were  temporarily  set  to  5.02988  km  and  0.015269  km/sec, 
respectively,  to  correspond  with  the  final  values  for  a  Clohessey-Wiltshire  rendezvous  at 
368  minutes.  The  values  in  the  control  weighting  matrix  were  then  adjusted  to  achieve 
the  Linear  Quadratic  Regulator  rendezvous  as  close  to  368  minutes  as  possible.  This 
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proved  difficult  because  of  the  sensitivity  of  the  controller.  Slightly  lower  control 
weighting  resulted  in  the  microsatellite  circling  the  target  satellite  more  times  as  it 
converged,  while  slightly  higher  control  weighting  produced  a  more  direct  flight  route. 
However,  the  difference  in  AH  between  the  more  and  less  direct  flight  paths  was 
minimal,  less  than  5  percent.  The  result  achieved  closest  to  368  minutes  was  321  minutes 
for  rendezvous,  which  required  0.098694  km/sec  of  AV  compared  to  0.016539  km/sec 
for  the  initial  AV  in  a  similar  Clohessey-Wiltshire  rendezvous. 

The  Linear  Quadratic  Regulator  rendezvous  criteria  were  returned  to  a  relative 
distance  of  1  meter  or  less  and  a  relative  velocity  of  1  cm/sec  or  less,  and  the  J2 
perturbation  was  added  producing  the  results  in  Table  8.  The  values  for  the  state 
weighting  matrix  and  the  control  weighting  matrix  were  1  and  1  e  1 3,  respectively,  and  the 
rendezvous  is  depicted  in  Figure  24. 


LQR  Rendezvous  -  J2  Perturbation 

AVrou,i  (km/sec) 

0.308920 

Rendezvous  Time  (min) 

658 

Table  8  LQR  Rendezvous  With  J2  Perturbation 
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Figure  24  Relative  Distance  During  LQR  Rendezvous  With  J2  Perturbation 

The  out  of  plane  effect  of  the  J2  perturbation  is  seen  to  be  the  most  time  consuming  for 
the  controller  to  overcome.  Countering  the  J2  perturbation  effect  only  increases  the  AV 
by  2.1  m/sec  over  the  non-perturbed  result.  Therefore,  the  values  corresponding  to  dz  in 
the  state  weighting  matrix  are  henceforth  multiplied  by  100  to  speed  transition  to  the  final 
state.  The  new  performance  values  are  shown  in  Table  9  and  the  effect  on  rendezvous 
dynamics  is  pictured  in  Figure  25. 


LQR  Rendezvous  -  J2  Perturbation 

mai  (km/sec) 

0.308844 

Rendezvous  Time  (min) 

509 

Table  9  LQR  Rendezvous  With  J2  Perturbation  and  Increased  dz  State  Weighting 
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Time  (minutes) 


Figure  25  Relative  Distance  During  LQR  Rendezvous  With  J?  Perturbation 
and  Increased  Sz  State  Weighting 

The  increased  weighting  on  the  5z  values  in  the  state  weighting  matrix  actually 
decreases  the  overall  AV  required  for  the  maneuver  slightly  by  decreasing  the  time  for 
rendezvous  149  minutes. 

Finally,  the  drag  perturbation  is  added.  The  additional  burden  proves  to  be  too  great 
for  the  controller  to  overcome  with  the  previous  control  weighting  matrix  values. 
Therefore,  control  emphasis  was  increased  by  lowering  the  values  in  the  control 
weighting  matrix  by  half,  from  lel3to5el2.  The  increase  in  control  emphasis  is 
sufficient  to  achieve  rendezvous,  and  results  in  a  AV  increase  of  more  than  75  m/sec 
over  the  previous  value,  while  the  time  to  achieve  rendezvous  is  reduced  by  nearly  2 
hours.  The  performance  values  are  found  in  Table  10  and  the  position  of  the 
microsatellite  relative  to  the  target  satellite  is  captured  in  the  Sr,  roS6  plane  in  Figure  26. 
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LQR  Rendezvous  -  J2  and  Drag  Perturbations 

mai  (km/sec) 

0.383119 

Rendezvous  Time  (min) 

384 

Table  10  LQR  Rendezvous  with  J2  and  Drag  Perturbations 


Figure  26  LQR  Rendezvous  with/?  and  Drag  Perturbations  in  Sr,roS6  plane 


Hybrid  Controller 

The  final  controller  considered  is  a  hybrid  controller.  The  Clohessey-Wiltshire 
controller  has  been  shown  to  be  useful  for  reducing  large  distances  with  low  fuel  usage, 
but  it  lacks  sufficient  robustness  to  complete  the  maneuver.  The  Linear  Quadratic 
Regulator  controller  has  been  shown  to  be  quite  robust,  but  it  uses  a  great  deal  of  fuel  in 
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reducing  large  distances.  The  logical  solution  is  to  combine  the  two  controllers  to  exploit 
their  strengths  and  minimizes  their  weaknesses. 

The  hybrid  controller  uses  the  Clohessey-Wiltshire  controller  to  calculate  the  initial 
impulsive  AV.  Then  the  target  is  pursued  until  the  specified  rendezvous  time.  Instead  of 
applying  a  second  AV  to  null  out  the  relative  velocity,  the  position  and  velocity  of  the 
microsatellite  relative  to  the  target  satellite  at  rendezvous  time  are  fed  directly  into  the 
Linear  Quadratic  Regulator  controller  as  starting  values.  The  LQR  controller  then 
completes  the  rendezvous  to  the  specified  relative  position  and  velocity  criteria.  Since 
the  LQR  controller  is  only  used  to  reduce  a  small  distance  at  the  end,  the  control 
emphasis  can  be  increased  without  a  significant  penalty  in  AV .  For  the  test  cases  in  this 
section,  the  control  weighting  matrix  values  are  set  to  lel2. 

As  with  the  previous  controllers,  the  first  test  case  of  the  hybrid  controller  involves  a 
rendezvous  with  no  perturbations.  The  performance  values  are  captured  in  Table  1 1  and 
the  rendezvous  is  depicted  in  Figure  27. 


Hybrid  Rendezvous  -  No  Perturbations 

AVCW  (km/sec) 

0.016539 

^ lqr  (km/sec) 

0.026642 

mat  (km/sec) 

0.043181 

Rendezvous  Time  (minutes) 

574 

Table  1 1  Hybrid  Rendezvous  With  No  Perturbations 
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Figure  27  Hybrid  Rendezvous  With  No  Perturbations 


Next,  the  J?  perturbation  is  added  to  the  scenario.  The  resultant  performance  values 
are  listed  in  Table  12  and  the  maneuver  is  illustrated  in  Figure  28.  Of  particular  interest 
is  the  dz  graph  which  shows  divergence  during  the  Clohessey-Wiltshire  portion  of  the 
rendezvous,  then  quickly  gets  damped  out  once  the  Linear  Quadratic  Regulator  controller 
takes  over. 


Hybrid  Rendezvous  -  J2  Perturbation 

AVCW  (km/sec) 

0.016578 

(km/sec) 

0.029722 

^ Total  (km/sec) 

0.046301 

Rendezvous  Time  (minutes) 

573 

Table  12  Hybrid  Rendezvous  With  J?  Perturbation 
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Figure  28  Hybrid  Rendezvous  With  J2  Perturbation 


Drag  is  added  for  the  final  test  case.  The  performance  values  are  listed  in  Table  13 
and  the  motion  of  the  microsatellite  relative  to  the  target  is  plotted  in  Figure  29. 


Hybrid  Rendezvous  -  J2  and  Drag  Perturbations 

A.VCW  (km/sec) 

0.016578 

^Vlqr  (km/sec) 

0.032306 

(km/sec) 

0.048885 

Rendezvous  Time  (minutes) 

590 

Table  13  Hybrid  Rendezvous  With  J2  and  Drag  Perturbations 
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V.  Conclusions  and  Recommendations 


The  important  values  from  the  test  cases  for  all  three  controllers  studied  are  compiled 
in  Table  14. 


|  Controller\Perturbations 

No  Perturbations 

.7?  Perturbation 

J2  and  Drag 

CW 

^ Total  (km/sec) 

0.031809 

0.032496 

0.035748 

Distance  (km) 

5.029880 

3.518942 

3.321252 

Time  (min) 

368 

368 

368 

LQR 

SVTl)lal  (km/sec) 

0.306821 

0.308844 

0.383119 

Distance  (km) 

<0.001 

<0.001 

<0.001 

Time  (min) 

510 

509 

384 

Hybrid 

AVromi  (km/sec) 

0.043181 

0.046301 

0.048885 

Distance  (km) 

<0.001 

<0.001 

<0.001 

Time  (min) 

574 

573 

590 

Table  14  Compilation  of  Test  Case  Results  For  All  Three  Controllers 


Conclusions 

It  can  be  readily  seen  that  adding  perturbation  effects  increases  the  needed  AV  for  all 
three  controllers.  It  is  also  clear  that  the  Clohessey-Wiltshire  controller  uses  the  least 
AV ,  but  does  not  achieve  the  desired  rendezvous  distance.  Accomplishing  multiple  CW 
maneuvers  is  not  effective  either,  especially  if  the  Sr  offset  is  non-zero.  This  can  be  seen 
in  Appendix  A,  where  other  examples  of  CW  results  over  a  range  of  Sr  and  roS0  values 
can  be  found. 
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The  Linear  Quadratic  Regulator  controller  achieves  rendezvous  to  the  desired  relative 
distance  and  velocity,  but  spends  a  great  deal  of  fuel  in  the  process.  Attempts  to  decrease 
fuel  usage  by  adjusting  the  state  weighting  or  control  weighting  matrices  lengthen  the 
amount  of  time  required  to  achieve  rendezvous,  and  if  taken  too  far,  could  result  in 
insufficient  control  usage  to  even  complete  the  rendezvous. 

The  Hybrid  controller  is  clearly  the  best  solution.  It  achieves  rendezvous  to  the 
specified  relative  distance  and  velocity  while  only  increasing  AV  by  approximately  50% 
over  the  corresponding  values  for  the  Clohessey-Wiltshire  controller. 

The  results  of  this  study  suggest  that  the  control  portion  of  the  microsatellite 
rendezvous  is  feasible  with  reasonable  fuel  usage. 


Recommendations 

In  the  test  cases  used  for  this  study,  the  microsatellite  began  rendezvous  from  1000 
kilometers  in-track  ( roS0 )  behind  the  target  satellite.  Further  study  should  include 
varying  the  starting  location  in  all  three  directions,  in-track,  radial  ( Sr  )  and  out  of  plane 
( Sz  ).  Some  examples  of  varying  starting  position  in  the  in-track  and  radial  directions  are 
found  in  Appendix  A. 

This  study  focused  on  satellites  in  Low  Earth  Orbit.  Future  study  should  expand  to 
include  higher  orbits.  Additionally,  this  study  included  a  small  eccentricity  in  the  target 
orbit.  Future  investigation  should  include  increasing  the  target  orbit’s  eccentricity  to 
explore  the  performance  envelope  of  the  controllers. 
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The  Clohessey-Wiltshire  equations  need  to  be  further  developed  to  account  for 
inclination  differences  between  the  microsatellite  and  the  target  satellite.  Development  of 
the  differential  inclination  terms  is  discussed  by  Swank,  et  al.  (8) 

The  greatest  potential  for  improvement  in  controller  perfonnance  is  incorporation  of 
gain  scheduling  in  the  Linear  Quadratic  Regulator.  The  control  usage  should  begin  low 
to  avoid  drastic  maneuvers  which  can  be  costly  in  fuel,  and  then  increase  as  the 
microsatellite  approaches  the  target,  to  complete  the  rendezvous.  To  implement  gain 
scheduling,  the  optimal  gain  matrix  could  be  recalculated  in  each  iteration  of  the 
propagation  algorithm.  The  control  weighting  could  be  set  proportional  to  the  remaining 
distance  between  the  microsatellite  and  the  target  satellite. 

This  study  was  conducted  as  though  perfect  knowledge  of  the  microsatellite’s  and  the 
target  satellite’s  position  and  velocity  were  known  at  all  times.  Future  study  should 
incorporate  realistic  uncertainties  in  these  values,  and  performance  of  the  controllers 
should  be  characterized  statistically.  Incorporation  of  a  sequential  filter  seems  to  be  the 
logical  next  step. 

Finally,  a  systems  engineering  study  should  be  conducted  to  determine  the  feasibility 
of  implementing  the  control  solutions  investigated  in  this  thesis.  The  propulsion  system, 
processor,  and  onboard  sensor  capability  should  be  key  areas  of  focus. 
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Appendix  A  Example  Controller  Results 


This  appendix  contains  tables  showing  controller  solution  results  over  a  range  of 
starting  offsets.  The  in-track  offset  ( ro86 )  is  varied  from  -1000  kilometers  to  +1000 

kilometers  while  the  radial  offset  ( Sr  )  is  varied  from  -1  kilometer  to  +1  kilometer.  No 
offsets  in  the  out  of  plane  direction  ( dz  )  are  included. 

Tables  15-17  are  the  results  of  using  the  Clohessey-Wiltshire  controller.  Table  15 
was  generated  without  perturbations;  Table  16  includes  the  J2  perturbation  and  Table  17 
adds  the  drag  perturbation.  All  Clohessey-Wiltshire  solutions  calculated  in  these  tables 
use  a  rendezvous  time  of  368  minutes,  though  it  should  be  clear  from  previous  discussion 
that  the  optimal  rendezvous  time  will  vary  with  the  relative  geometry  of  the  microsatellite 
and  target  satellite  to  the  circular  reference  frame.  The  top  value  in  each  cell  is  the  total 
AV  (km/sec)  needed  for  the  maneuver,  and  the  lower  is  the  relative  distance  (km)  at  the 
specified  rendezvous  time. 

Tables  18-20  are  the  results  of  using  the  Linear  Quadratic  Regulator  controller. 

Table  18  was  generated  without  perturbations;  Table  19  includes  the  J2  perturbation  and 
Table  20  adds  the  drag  perturbation.  All  LQR  solutions  were  calculated  with  state 
weighting  matrix  values  of  1.  Table  18  and  Table  19  were  calculated  using  control 
weighting  matrix  values  of  lei 3  while  the  control  weighting  matrix  value  used  for 
Table  20  is  le  12.  The  top  value  in  each  cell  represents  the  total  AV  (km/sec)  for  the 
maneuver.  The  lower  number  is  the  time  that  it  takes  in  integer  minutes  for  the 
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microsatellite  to  rendezvous  to  within  the  specified  criteria  for  relative  distance  (<  lm) 
and  relative  velocity  (<  lcm/sec)  of  the  target  satellite. 

Tables  21-23  are  the  results  of  using  the  Hybrid  controller.  Table  21  was  generated 
without  perturbations;  Table  22  includes  the  J2  perturbation  and  Table  23  adds  the  drag 
perturbation.  For  all  Hybrid  controller  solutions  in  these  tables,  the  CW  portion  was 
calculated  using  a  rendezvous  time  of  368  minutes,  and  the  LQR  solutions  used  state 
weighting  matrix  values  of  1  and  control  weighting  matrix  values  of  le  12.  The  top  value 
in  each  cell  represents  the  total  AV  (km/sec)  for  the  maneuver.  The  lower  number  is  the 
time  that  it  takes  in  integer  minutes  for  the  microsatellite  to  rendezvous  to  within  the 
specified  criteria  for  relative  distance  (<  lm)  and  relative  velocity  (<  lcm/sec)  of  the 
target  satellite. 
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-1000 

-100 

-10 

-1 

0.03445 

0.01857 

0.01946 

0.01958 

5.58052 

0.84143 

1.15968 

1.19841 

0.03159 

0.00318 

0.00182 

0.00191 

5.05209 

0.16314 

0.06485 

0.09129 

0.03180 

0.00332 

3.35E-4 

3.35E-5 

5.02988 

0.24071 

0.03131 

0.00320 

0.03212 

0.00439 

0.00206 

0.00193 

5.01508 

0.31387 

0.11700 

0.09181 

0.03873 

0.02019 

0.01896 

0.01885 

5.21682 

0.74610 

0.68924 

0.68031 

0.01960 


1 .20280 


0.00192 


0.09439 


0.00191 


0.08902 


0.01884 


0.67930 


1 

10 

100 

1000 

0.01961 

0.01974 

0.02127 

0.05084 

1.20719 

1 .24745 

1.71914 

13.8511 

0.00193 

0.00207 

0.00451 

0.03738 

0.09752 

0.12681 

0.50169 

11.7393 

3.35E-5 

3.35E-4 

0.00337 

0.03619 

0.00322 

0.03293 

0.40283 

11.5476 

0.00190 

0.00181 

0.00313 

0.03507 

0.08623 

0.06210 

0.31697 

11.3650 

0.01882 

0.01872 

0.01783 

0.02956 

0.67829 

0.66904 

0.59703 

10.1483 

Table  15  CW  Solutions  Without  Perturbations 


-1000 

-100 

-10 

-1 

1 

10 

100 

1000 

1 

0.03292 

0.01764 

0.01979 

0.02005 

0.02008 

0.02011 

0.02037 

0.02339 

0.07507 

4.27274 

0.23973 

1.03635 

1.17372 

1.18919 

1.20469 

1.34594 

2.92660 

35.1488 

D 

0.03202 

0.00340 

0.00171 

0.00194 

0.00197 

0.00200 

0.00227 

0.00594 

0.05954 

3.55563 

1.21481 

0.06536 

0.07447 

0.08863 

0.10301 

0.23663 

1.74987 

33.1284 

0 

0.03249 

0.00420 

4.33E-4 

4.34E-5 

4.35E-5 

4.36E-4 

0.00448 

0.05802 

3.51894 

1.29209 

0.14557 

0.01472 

0.01475 

0.14921 

1.65601 

32.9425 

ra 

0.03307 

0.00560 

0.00226 

0.00199 

0.00196 

0.00193 

0.00170 

0.00357 

0.05655 

3.49129 

1.36253 

0.22425 

0.09607 

0.08204 

0.06822 

0.07474 

1.56994 

32.7643 

-i 

0.04170 

0.02206 

0.01956 

0.01931 

0.01928 

0.01925 

0.01901 

0.01663 

0.04608 

3.65082 

1.68858 

0.65572 

0.54851 

0.53689 

0.52536 

0.42699 

1.17010 

31.5111 

Table  16  CW  Solutions  With  J2  Perturbations 


-1000 

-100 

-10 

-1 

1 

10 

100 

1000 

1 

0.03180 

0.01396 

0.01628 

0.01655 

0.01659 

0.01662 

0.01691 

0.02017 

0.07335 

2.69499 

3.47650 

2.11697 

1.97274 

1.95656 

1.94036 

1.79366 

0.67160 

33.5371 

D 

0.03490 

0.00646 

0.00357 

0.00340 

0.00339 

0.00337 

0.00324 

0.00399 

0.05787 

3.27023 

4.50285 

3.15144 

3.00127 

2.98441 

2.96753 

2.81408 

1.14845 

31.5088 

0 

0.03574 

0.00757 

0.00389 

3.53E-3 

3.49E-3 

0.00343 

0.00366 

0.05635 

3.32125 

4.58521 

3.23852 

3.08872 

3.05507 

2.90196 

1.23276 

31.3205 

a 

0.03666 

0.00911 

0.00577 

0.00549 

0.00546 

0.00543 

0.00517 

0.00418 

0.05489 

3.36905 

4.66120 

3.31944 

3.17007 

3.15330 

3.13650 

2.98379 

1.31361 

31.1397 

-i 

0.04745 

0.02585 

0.02314 

0.02288 

0.02285 

0.02282 

0.02255 

0.01993 

0.04448 

3.65081 

5.05811 

3.76722 

3.62303 

3.60684 

3.59062 

3.44309 

1.81961 

29.8511 

Table  17  CW  Solutions  With  J?  and  Drag  Perturbations 
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-1000  -100  -10  -1  0  1  10  100  1000 

0.30791  0.03211  0.00428  0.00156  0.00128  0.00102  0.00210  0.03029  0.33944 


510 

425 

347 

286 

282 

280 

310 

424 

506 

n 

0.30693 

0.03110 

0.00322 

4.28E-4 

1.28E-4 

2.10E-4 

0.00300 

0.03129 

0.34050 

510 

425 

318 

241 

199 

228 

317 

424 

506 

0 

0.30682 

0.03099 

0.00311 

3.1  IE-4 

3.1  IE-4 

0.00312 

0.03140 

0.34062 

510 

424 

317 

235 

235 

317 

424 

506 

ra 

0.30671 

0.03088 

0.00300 

2.10E-4 

1.28E-4 

4.28E-4 

0.00323 

0.03152 

0.34074 

510 

424 

317 

228 

199 

241 

318 

424 

506 

-i 

0.30572 

0.02988 

0.00210 

0.00102 

0.00128 

0.00156 

0.00428 

0.03253 

0.34180 

510 

424 

310 

280 

282 

286 

347 

425 

506 

Table  18  LQR  Solutions  Without  Perturbations 


-1000 

-100 

-10 

-1 

1 

10 
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1000 
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0.30994 

0.03199 

0.00426 

0.00156 

0.00128 
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0.00210 

0.03015 

0.33706 
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280 
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509 
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317 
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423 
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0.03088 

0.00310 
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0.03126 
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423 
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ra 

0.30873 

0.03077 

0.00299 

2.09E-4 

1.28E-4 

4.26E-4 

0.00321 

0.03137 

0.33835 

509 

423 

316 

227 

199 

240 

317 

423 

505 

-i 

0.30775 

0.02978 

0.00209 

0.00102 

0.00128 

0.00156 

0.00426 

0.03238 

0.33940 

509 

423 
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280 
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346 

424 
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Table  19  LQR  Solutions  With  J?  Perturbation 
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0.70570 

0.07758 
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0.00269 

0.00202 
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0.07695 
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0 

0.70399 

0.07580 
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0.00767 

0.07713 

0.84045 
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181 
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H 
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5.95E-4 

2.12E-4 

9.60E-4 

0.00784 

0.07731 

0.84064 

301 

239 

181 

120 

99 

146 

207 

218 

275 

-i 

0.70228 

0.07403 

0.00592 

0.00133 

0.00198 

0.00269 

0.00948 

0.07892 

0.84231 

301 

239 

179 

155 

154 

155 

208 

218 

275 

Table  20  LQR  Solutions  With  J2  and  Drag  Perturbations 


66 


-1000  -100  -10  -1  0  1  10  100  1000 

0.03834  0.02072  0.02210  0.02227  0.02228  0.02230  0.02247  0.02442  0.05463 


577 

563 

564 

564 

564 

564 

564 

564 

567 

n 

0.04214 

0.00364 

0.00204 

0.00219 

0.00221 

0.00223 

0.00243 

0.00555 

0.04207 

574 

520 

506 

507 

507 

507 

507 

529 

549 

0 

0.04318 

0.00415 

4.12E-4 

2.55E-5 

2.55E-5 

4.1  IE-4 

0.00413 

0.04097 

574 

513 

452 

374 

374 

451 

512 

549 

ra 

0.04431 

0.00558 

0.00243 

0.00223 

0.00221 

0.00219 

0.00204 

0.00361 

0.03995 

574 

531 

507 

507 

507 

507 

506 

519 

549 

-i 

0.05785 

0.02438 

0.02220 

0.02200 

0.02198 

0.02196 

0.02177 

0.02017 

0.03552 

575 

565 

564 

564 

564 

564 

564 

563 

560 

Table  21  Hybrid  Solutions  Without  Perturbations 


-1000 

-100 

-10 

-1 

1 

10 

100 

1000 

1 

0.03925 

0.02071 

0.02285 

0.02310 

0.02312 

0.02315 

0.02340 

0.02621 

0.07039 

575 

563 

564 

564 

564 

564 

564 

565 

579 

0.1 

0.04505 

0.00388 

0.00204 

0.00226 

0.00229 

0.00232 

0.00261 

0.00657 

0.05644 

573 

457 

506 

507 

507 

508 

509 

531 

561 

0 

0.04630 

0.00485 

4.85E-4 

4.80E-5 

4.80E-5 

4.86E-4 

0.00494 

0.05513 

573 

506 

440 

425 

425 

440 

503 

560 

ra 

0.04763 

0.00649 

0.00261 

0.00232 

0.00229 

0.00226 

0.00203 

0.00396 

0.05388 

573 

532 

509 

508 

507 

507 

506 

497 

559 

-i 

0.06262 

0.02605 

0.02310 

0.02283 

0.02280 

0.02277 

0.02250 

0.02007 

0.04596 

586 

565 

564 

564 

564 

564 

564 

563 

553 

Table  22  Hybrid  Solutions  With  J?  Perturbation 


-1000 

-100 

-10 

-1 

1 

10 

100 

1000 

1 

0.03920 

0.01904 

0.02125 

0.02150 

0.02153 

0.02156 

0.02182 

0.02465 

0.06925 

549 

562 

562 

562 

562 

562 

562 

563 

561 

ra 

0.04747 

0.00611 

0.00289 

0.00276 

0.00275 

0.00274 

0.00270 

0.00549 

0.05531 

591 

520 

517 

516 

516 

516 

516 

514 

556 

0 

0.04888 

0.00725 

0.00309 

2.71E-3 

2.66E-3 

0.00260 

0.00432 

0.05400 

590 

524 

521 

520 

520 

520 

517 

556 

ra 

0.05036 

0.00892 

0.00504 

0.00472 

0.00469 

0.00466 

0.00438 

0.00416 

0.05276 

590 

527 

525 

524 

524 

524 

524 

521 

555 

-i 

0.06608 

0.02817 

0.02499 

0.02469 

0.02466 

0.02462 

0.02433 

0.02170 

0.04489 

590 

585 

584 

584 

584 

584 

584 

584 

546 

Table  23  Hybrid  Solutions  With  J2  and  Drag  Perturbations 
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Appendix  B  Linear  Quadratic  Regulator  Propagation  Algorithm 
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Appendix  C  MATLAB  Code 


MAIN 


% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

%■ 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

%= 


THESIS  -  MAIN  CODE 
AS  OF:  04  MAR  03 
Troy  Tschirhart 

Orbital  Rendezvous  With  a  Non-Cooperative  Target 


This  program  uses  the  following  function  files  which  must  be  on  the  current  path: 


atmosphere.m 

Calclnit.m 

CWRend 

cwvijk 

DoPlots 

LQRRend 

propagate,  m 

posvel.m 

ijk2pqw.m 

pqw2ijk.m 

rtz2pqw.m 

rv2coe.m 

coe2rv.m 


calculate  atmospheric  density  at  the  given  altitude 
calculate  the  initial  conditions  for  the  run 
accomplish  clohessey-wiltshire  rendezvous  maneuver 
calculate  the  clohessey-wiltshire  velocity  in  inertial  frame 
plot  the  results 

accomplish  lqr  rendezvous  manuever 
propagator 

set  up  the  differential  equation  for  the  propagator 
transform  r,v  from  ijk  frame  to  pqw  frame 
transform  r,v  from  pqw  frame  to  ijk  frame 
transform  r,v  from  rtz  frame  to  pqw  frame 
calculate  coe  for  the  given  r,v 
calculate  r,v  for  the  given  coe 
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%============================= 

% 

%  Clear  Variables  and  Set  Format  Options 

% 

%============================= 


clear 

format  long  g 
format  compact 


%================================================= 

% 

%  Print  a  banner  to  separate  results 

%  Start  the  timer  (used  at  the  end  to  determine  how  long  the  run  took) 

% 

%================================================= 

((================================================== 


0 


tic 


%======— ========= 

% 

%  Set  Selectable  Variable  Values 

% 

%======_========== 


% - 

%  The  target's  actual  initial  COEs 
% - 


coe_tgt_act(l) 

coe_tgt_act(2) 

coe_tgt_act(3) 

coe_tgt_act(4) 

coe_tgt_act(5) 

coe_tgt_act(6) 

coe_tgt_act(7) 


6.7728889 12204840e+003; 

%  km 

a 

9. 8877 1 35498259 13e-004; 

%  dimensionless 

e 

0.79736386485827; 

%  radians 

nu 

0.90757990078380; 

%  radians 

i 

1.51843760980691; 

%  radians 

capomega 

5.59054044657763; 

%  radians 

small  omega 

0.0; 

%  seconds 

time  since  perigee 
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% - 

%  Set  initial  micro  offset  from  the  target 
% - 

dist  =  -1000;  %  kilometers  arclength  (ro*delta_theta) 

delr  =  0;  %  kilometers  (delta_r) 

delz  =  0;  %  kilometers  (delta  z) 

% - 

%  Set  the  acceptable  relative  distance  and  velocity  for  a  successful  rendezvous 
% - 

catchdis  =  0.00 1 ;  %  kilometers 

catchvel  =  0.00001;  %  kilometers/second 


% - 

%  Set  controller  options  (Notes:  "1"  =  option  selected;  "0"  =  option  not  selected; 

%  Only  one  controller  option  below  should  be  selected  for  each  run) 

% - 

CW  =1;  %  Use  the  Clohessey-Wiltshire  controller  only 

CWLR  =  0;  %  Use  the  Clohessey-Wiltshire  -  LQR  hybrid  controller 

LR  =  0;  %  Use  the  LQR  controller  only 


% - 

%  Set  step  size  and  CW  rendezvous  time 
% - 


timestep  =  60; 

%  seconds 

rendtime  =  368; 

%  Integer  number  of  timesteps 
%  Examples:  For  a  rendezvous  time  of  184  minutes 
%  184  if  timestep  =  60;  1 1040  [=184*60]  if  timestep  =  1 
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% - 

%  Specify  values  for  the  state  weighting  matrix,  Q 
%  and  the  control  weighting  matrix,  R 

%  Note:  Q_mag  increases  =>  faster  movement  from  initial  to  desired  states 
%  R_mag  increases  =>  lower  control  usage 
% - 


Q_mag  =  1; 
Rmag  =  lel3; 


% - 

%  Set  perturbation  (J2)  option  (Note: "  1"  =  option  selected;  "0"  =  option  not  selected) 
% - 

pert  =  0;  %  Include  perturbation  terms 


% - 

%  Set  drag  options  and  values  (Note:  "1"  =  option  selected;  "0"  =  option  not  selected) 
% - 

dragtgt  =  0;  %  Include  drag  in  target's  propagations 

cd_tgt  =  2.2;  %  Drag  coefficient  of  the  target 

a_tgt  =  3.5*  1 .2/(1000A2);  %  Area  of  the  target  (kmA2) 

m_tgt  =  725;  %  Mass  of  the  target 

cdamtgt  =  dragtgt  *  (cd_tgt  *  a_tgt)  /  m_tgt;  %  Calculate  the  target's  cdam  value 

% - 


dragmic  =  0; 
cdmic  =  3; 
amic  =  1.5/(1000A2); 
mmic  =100; 
cdammic  =  dragmic  *  (cd  mic  *  a  mic)  /  m 


%  Include  drag  in  micro's  propagations 
%  Drag  coefficient  of  the  micro 
%  Area  of  the  micro  (kmA2) 

%  Mass  of  the  micro 
l  %  Calculate  the  micro's  cdam  value 
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% - 

%  Set  plot  options  (Note:  "1"  =  option  selected;  "0"  =  option  not  selected) 

% - 

prdijk  =  1 ;  %  Plot  relative  distance  in  the  inertial  (ijk)  frame 

prdrtz  =  1 ;  %  Plot  relative  distance  in  the  relative  (rtz)  frame 

prdroto  =1;  %  Plot  relative  distance  in  the  relative  plane  (delta  r,  ro*delta_theta) 


%======—====—======— ============= 

% 

%  Initialize  variable  values 

% 

%========================================== 

points  =  rendtime  +  1 ;  %  set  the  number  of  points  for  propagation 

%  (note;  point  1  is  really  0) 

deltavaccum  =  0;  %  initialize  delta-V  to  zero 


%===_======_=====__= 

% 

%  Calculate  Initial  Values 

% 

%  1.  Target's  initial  position  and  velocity 
%  2.  Micro's  initial  position  and  velocity 

% 

%======—=====—======= 

Calclnit 
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%============================================================= 

% 

%  Accomplish  Clohessy-Wiltshire  Rendezvous 

% 

%  1 .  Determine  the  micro's  position  relative  to  the  circular  reference  at  the  start  time 
%  2.  Determine  the  target's  position  relative  to  the  circular  reference  at  CW  rend  time 
%  3.  Solve  the  CW  equations 

%  4.  Transform  the  delta-V  from  the  relative  frame  to  the  inertial  frame 
%  5.  Apply  the  delta-V  to  the  micro 
%  6.  Fly  out  the  micro  and  the  target  to  rendezvous  time 
%  7.  Check  final  relative  distance  and  velocity  between  the  micro  and  the  target 

% 

%—— - - 


if  (CW  1)  |  (CWLR  ==  1) 

CWRend 
end 


%============================================================== 

% 

%  Accomplish  Linear  Quadratic  Regulator  Rendezvous 

% 

%  1 .  Starting  values:  position  and  velocity  of  the  target,  ref,  and  micro  at  CW  rend  time 
%  2.  Calculate  the  A  matrix,  set  up  the  B,  Q,  and  R  matrices 
%  3.  Solve  for  the  gain  matrix,  K 
%  4.  Iterate 

%  A.  Calculate  the  necessary  control  input 
%  B.  Propagate  one  step 

%  C.  Repeat  until  successful  rendezvous  is  achieved 

% 

%==—==—=—=—=—===—=—== 


if  LR  ==  1 
LQRRend 
elseif  CWLR  ==  1 

rmicroijktmp  =  r_micro_ijk(rend_time,:); 
vmicroijktmp  =  v_micro_ijk(rend_time,:); 
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rtgtacttmp  =  r_tgt_act_ijk(rend_time,:); 
v_tgt_act_tmp  =  v_tgt_act_ijk(rend_time,:); 
rrefijktmp  =  r_ref_ijk(rend_time,:); 
v_ref_ijk_tmp  =  v_ref_ijk(rend_time,:); 
LQRRend 
end 


%================ 

% 

%  Print  Output  Values 

% 

%================ 


deltar  =  delr 
rodeltatheta  =  dist 
deltaz  =  delz 

ifCW  ==  1 


pert 

dragtgt 

dragmic 

rendtime 

timestep 

delta_v_cwl 

delta_v_cw2 

deltavaccum 

dist_at_rend_time  =  r_mag_rel_cw 

elseif  CWLR  ==  1 

pert 

dragtgt 

dragmic 

rendtime 

pt 

rend_time+pt 

timestep 
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Qmag 

Rmag 

delta_v_cwl 

delta_v_lqr 

deltavaccum 

dist_at_rend_time  =  r_mag_rel_cw 
finalvel  =  vel_now 
finaldist  =  distnow 

elseif  LR  ==  1 


pert 

dragtgt 

dragmic 

pt 

timestep 

Qmag 

Rmag 

delta_v_lqr 

delta_v_accum 

finalvel  =  vel_now 

finaldist  =  distnow 

end 


%============— 

% 

%  Draw  Desired  Plots 

% 

%============— 

Do  Plots 


%============= 

% 

%  End  of  Program 

% 

%— ========== 

run  time  =  toe 
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